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ABSTRACT 

The current dynamical structure of the Kuiper belt was shaped by the orbital evolution of the 
giant planets, especially Neptune, during the era following planet formation, when the giant plan- 
ets may have undergone planet-planet scattering and/or planetesimal-driven migration. Numerical 
simulations of this process, while reproducing many properties of the belt, fail to generate the high 
inclinations and eccentricities observed for some objects while maintaining the observed dynamically 
"cold" population. We present the first of a three-part parameter study of how different dynamical 
histories of Neptune sculpt the planetesimal disk. Here we identify which dynamical histories allow an 
in situ planetesimal disk to remain dynamically cold, becoming today's cold Kuiper belt population. 
We find that if Neptune undergoes a period of elevated eccentricity and/or inclination, it secularly 
excites the eccentricities and inclinations of the planetesimal disk. We demonstrate that there are 
several well-defined regimes for this secular excitation, depending on the relative timescales of Nep- 
tune's migration, the damping of Neptune's orbital inclination and/or eccentricity, and the secular 
evolution of the planetesimals. We model this secular excitation analytically in each regime, allowing 
for a thorough exploration of parameter space. Neptune's eccentricity and inclination can remain 
high for a limited amount of time without disrupting the cold classical belt. In the regime of slow 
damping and slow migration, if Neptune is located (for example) at 20 AU, then its eccentricity must 
stay below 0.18 and its inclination below 6°. 

Subject headings: Kuiper Belt, planets and satellites: Neptune, solar system: general 



1. INTRODUCTION 

The solar system is often used as a case study for 
the formation of planetary systems from proto-planetary 
disks. The current configuration of Kuiper belt objects 
(KBOs) provides a map for how the dynamical evolu- 
tion of the giant planets in our solar system sculpted the 
disk of planetesimals. Therefore it is possible to use the 
orbital properties of the Kuiper belt to constrain how 
the orbits of the giant planets evolved in the early so- 
lar system, particularly for Neptune, the primary sculp- 
tor of the Kuiper belt. Models employing N-body in- 
tegrations to trace the effects of the giant planets' mi- 
gration and orbital eccentricity evolution on the plan- 
etesimal disk have enjoyed substantial success in repro- 
ducing t he dynamical populations of KBOs observed to- 
day (e. g . lMalhotralll993L 1199ft lHahn fc Maihotral 119991: 
Gomesl 120031: lHahn fc Malhotral 1200ft iLevison et al.l 
2008HMorbidelli et al.ll2008t ). These populations include 
objects near orbital resonance with Neptune, objects 
scattering off Neptu ne, and "classical" obj ects decoupled 
from Neptune. (See lGladman et al.l l2008. for definitions 
of the dynamical classes.) Yet substantial discrepancies 



still exist between the simulations and observations, par- 
ticularly for the classical population. The bias-corrected 
inclina tion distribution of observed classical KBOs is bi- 
moda l (IBrownll200lHGulbis et al.1l2010l:IVolk fc Malhotral 
1201 If) . The low- inclination, dynamically "cold" compo- 
nent and the high-inclination, dynamically "hot" com- 
ponent al so have distinct physical properties, includ- 
ing colors dTeeler fc Romani shin 200 0tlTruiillo fc Brownl 
[20031: iPeixinho et al.ll200 8fl. sizes (ILevison fc Sternl l200H; 
iFraser et al.ll2010D. albedos (iBrucker et al.ll2009D. and b i- 
nary fractions (jStephens fe Nollll2006HNoll et al.ll2008[) . 
To date no simulations have been able to produce both 
the high and low inclination classical objects while qual- 
itatively matching their observed eccentricity distribu- 
tion. 

Several theories of the dynamical history of the giant 
planets in our solar system have been proposed, inspired 
by the dynamical populations within the Kuiper belt. 
One of the most widely accepted models, the Nice Model, 
stems from a postul ated large scale instabi li ty in the early 
solar system (e.g. IThommes et al.l 119991 : iGomed 120031 : 
Morbidelli et al.l 12008ft . Inspired by the Nice Model, 



Levison et al. (2008) proposes a scenario in which Nep- 
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tune is scattered outward onto an eccentric orbit to near 
its current lo cation from a formati o n loc ation closer to 
the sun (e.g. IThommes et~all [r999. 2002). In this see- 
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nario, Neptune's eccentricity subsequently damps due to 
dynamical friction with a remnant disk of planetesimals. 
This model reproduces the observed resonant popula- 
tion, scattered population, and the hot classical popula- 
tion, but, we argue, does not satisfactorily match the low 
eccentricities of the observed population of dynamically 
cold objects in the classical region. Furthermore, it does 
not produce a sufficient number of high-inclinat ion ob- 
jects. Other incarnations of the Nice model (e.g. iGomesI 
120031: iMorbidelli et alJ l2008f > include an in situ popula- 
tion of cold objects, but, over the course of Neptune's 
evolution, these objects become excited to higher eccen- 
tricities. Whether or not the particular history described 
by the Nice Model is correct, the prevalence of large ec- 
centricities among extrasolar giant planets suggests that 
large-scale orbital excitation is common during the for- 
mation of planetary systems. 

In co ntrast to the uph eaval of the Nice Model, another 
model (|MalhotralfT993l) proposes a period of extensive, 
smooth migratio n of the giant planets . Planetesimal- 
driven migration ([Fernandez fe Ip| |l984) is likely impor- 
tant i n shaping t he architect u re of many planeta r y sys - 
tems. IMalhotral ([1993L 119951 ): lHahn fc Malhotral (1999) 
demonstrated that the migration of the giant planets 
in our solar system on flat, circular orbits would have 
perturbed the disk of planetesimals, scattering some to 
more eccentric orbits and capturing others into reso- 
nance. They also found that migration results in a large 
number of planetesimals in the location known today as 
the "Scattered Disk," which encompasses those objects 
that have had gravitational interactions with one or more 
of the giant planets and have high orbital eccentricities 
and inclinations. The key inconsistency between this 
model and current observations is that it cannot pro- 
duc e, without additional p roces ses such as stochasticity 
(e.glLevison fe Morbide"Iiill2003l : iMurrav-Clav fc Chiang! 
I2006f) . both the cold and hot classical populations from 
a single set of initial disk conditions or account for the 
differences in their physical properties. 

Because detailed N-body simulations are computation- 
ally expensive, previous works have investigated a limited 
number of solar system planetary histories. Currently, no 
comprehensive parameter study has been done exploring 
the evolution of Neptune's orbit given constraints from 
the sculpting of the Kuiper belt. In addition, although 
planet-pla net scattering often pro duces large mutual in- 
clinations (jChatteriee et al.ll2011lh as in observed in the 
Upsilon Andromedae svstem ()McArthur et al.l[2010l ). no 
serious, detailed treatment has included the possibility 
that Neptune underwent a period of high orbital inclina- 
tion, a conceivable outcome of the orbital instability in 
the early solar system. Here we consider a general model 
that can encompass the two detailed models described 
above - the Nice-Model-inspired scenario in which Nep- 
tune is scattered to a high eccentricity ([Levison et al.l 
120081 ) and the scenario of extensive migration o f Neptune 
on a low-eccentricity, low-inclination orbit (e.g.. IMalhotral 
1995) - as well as potential scenarios in which Neptune 
undergoes a period of high inclination. In this gener- 
alized model, Neptune undergoes some combination of 
migration and/or evolution of its eccentricity and/or in- 
clination. Here we take a step toward understanding the 
qualitative differences in the dynamics of the Kuiper belt 
generated by a wide range of planetary histories. In a 



series of three papers, we perform a parameter study of 
the effects on a disk of planetesimals of the migration and 
the eccentricity and inclination evolution of Neptune in 
the early solar system. No component of this parame- 
terization is new; rather, our goal is to comprehensively 
consider all possible parameters for Neptune's history. 

In this first paper, we introduce our method for com- 
putationally and analytically modeling Neptune's orbital 
evolution and its effects on the planetesimal population. 
Then we demonstrate several concepts that will allow us 
to thoroughly explore the parameter space of Neptune's 
dynamical history: 

• If Neptune undergoes a period of elevated orbital 
eccentricity, it will secularly excite the eccentrici- 
ties and inclinations of an in situ planetesimal pop- 
ulation. We can analytically model this secular ex- 
citation using a simple expression. 

• As Neptune's eccentricity and inclination damp, 
the planetesimals evolve to final eccentricities and 
inclinations that depend on Neptune's initial ec- 
centricity and inclination and damping timescales. 

• The effects of Neptune's eccentricity and inclina- 
tion evolution on the planetesimals can be treated 
separately to first order. 

• The migration rate sets Neptune's effective location 
for the secular evolution of the planetesimals. 

Having demonstrated these points, we can place robust 
constraints on Neptune's semi-major axis, eccentricity, 
and inclination during its late evolution (after any period 
of planet-planet scattering) , and on its migration, eccen- 
tricity damping, and inclination damping rates, identify- 
ing which parameters are consistent with maintaining the 
low eccentricities an d inclinations of the cold classi cals. 
In the second paper (|Dawson fc Murrav-Clavll2012| ), we 
place even stronger constraints by incorporating addi- 
tional effects, including the effects of other planets on 
Neptune's orbital evolution and the effects of proximity 
to mean-motion resonance with Neptune on the secular 
excitation of the planetesimals, as well as additional con- 
straints from the hot classical population. A third paper 
(Dawson and Murray-Clay 2012b, in prep) will focus on 
the inclinations of the classicals. 

In Section [2[ we present our parameterization of Nep- 
tune's orbital evolution. In Section [3l we discuss the ob- 
servational constraints that today's cold classical KBOs 
place on past sculpting of the planetesimal disk and de- 
scribe our computational and analytical models of the 
excitation of the planetesimal disk by an inclined and 
eccentric Neptune. We present the results of the param- 
eter study in Section 01 In Section [5j we present our 
conclusions and describe how the companion papers will 
expand upon this work. 

2. MODELING NEPTUNE'S ORBITAL EVOLUTION 

As a first step toward a comprehensive study of the 
impact of a wide range of dynamical histories of the 
outer solar system on the Kuiper belt, we thoroughly 
explore the parameter space of a general model for Nep- 
tune's dynamical history. This model encompasses spe- 
cific, previously-proposed solar system history models 
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(e.g. Malhotralll995l: Eevison et alJ ^OOS). In our param- 
eterization, Neptune is instantaneously scattered to an 
initial location with an initial eccentricity and inclina- 
tion. Subsequently, it undergoes planetesimal-driven mi- 
gration and its eccentricity and inclination are damped 
by dynamical friction from the disk of planetesimals. 
Resonant relaxation may also contribute significantly to 
the planet's eccentricity and inclination damping if the 
Toomre parameter of the planetesimal disk is large or 
if the random v elocities of the planetesimals are large 
(Tr cmainc 1998). The parameters of this model arc Nep- 
tune's "initial" (defined below) semi-major axis, eccen- 
tricity and inclination, the planet's migration rate, and 
the timescales for Neptune's eccentricity and inclination 
damping. In this model, we include the effects of only one 
planet (Neptune), an approach we justify briefly in Sec- 
tion !3.3l and more thoroughly in lDawson fc Murray-Clayl 
(|2012f ). 

2.1. Parameters 

We define Neptune's orbital evolution using the follow- 
ing parameters: 

• An "initial" semi-major axis, eccentricity, and in- 
clination. The initial semi-major axis is not nec- 
essarily the location where Neptune formed. We 
imagined that Neptune underwent a scattering, or 
series of scatterings, onto an inclined and eccentric 
orbit, which happened quickly enough to not dis- 
rupt the cold classical belt. We model the period 
after the scatterings end and Neptune's eccentricity 
and inclination begin to damp. 

• A migration rate, defined in Section [2.21 

• An eccentricity damping rate and inclination 
damping rate, also defined in Section l2~2l 

2.1.1. Values for the migration parameters 

Here we consider what range of parameters we should 
explore for Neptune's migration direction, distance, and 
timescale: 

• Migration direction: We know that Neptune's mi- 
gration will be, on average, outward because Nep- 
tune is exterior to a very massive planet (in this 
case Jupiter). As Neptune scatters planetesimals, 
Jupiter ejects them, leading to a net loss in an- 
gular momentum for Jupite r and gain in angular 
momentum for Neptu ne (e.g. lFernandez fe Ip|ll984t 
lMalhotralfT993lll995D . 

• Migration distance: If all the KBOs in orbital res- 
onance with Neptune were captured during migra- 
tion and from orbits with low eccentricity, Nep- 
tune needs to have migrated a distance of 7-10 AU 
l|Malhotralll99l fl99H IHahn fc Malhotrall2005t) to 
adiabatically raise their eccentricities to the val- 
ues observed. However, other mechanisms have 
been proposed for producing the population of 
resonant KBOs, in which case the 7-10 AU con- 
straint on the migra tion distance would not apply. 
iLevison et all (2008) argue that resonant objects 



were scattered from the inner disk into the clas- 
sical region, entered resonances widened by Nep- 
tune's high eccentricity, and were trapped when 
Neptu ne's eccentricity damped . In a companion 
paper ([Dawson fc Murrav-Clavll2012f) , we demon- 
strate that objects scattered into classical region 
secularly evolve very quickly near resonances, al- 
lowing them to reach stable, lower-eccentricity or- 
bits. Future observations of the binary fractions 
(and pos sibly colors) of resonant KBOs - p roper- 
ties that IMurrav-Clav fc Sc hlichting (2011) argue 
may encode the location of each KBO's formation - 
may distinguish between these mechanisms. There- 
fore we consider all "initial" semi-major axes for 
Neptune interior to its current location. We quote 
constraints for 20 AU and 30 AU as examples of a 
long and short migration distance, respectively. 

• Migration timescale: The distribution of libration 
angles of twotinos, KBOs in the 2:1 resonance, 
place a lower limi t of 1 Myr on Neptune's m i- 
gration timescale ([Murray- Clay fc Chiang! 120051) . 
(Upcoming unbiased surveys, e.g. PAN-STARRS, 
LSST, are needed to confirm the true distribu- 
tion of the twotinos.) An upper limit could be 
imposed by the stochasticity of the migration, 
which due to the finite sizes of planetesimals, 
limits the efficiency of keeping obj ects in reso- 
nance (|Murrav-Clav fc Chianel l2006f) . However, 
the amount of stochasticity depends on the size dis- 
tribution of the planetesimals, which is unknown 
and depends on the phy sics of their formation (see 
IChiang fc Youdinl 120101 and references therein). 
Given these uncertainties, we consider all migra- 
tion timescales greater than 0.3 Myr. 

2.2. Computational model 

Direct computational modeling of the effect of plan- 
etesimals on Neptune's orbit would be computation- 
ally expensive, so instead we apply fictitious forces 
(Appendix) to evolve Neptune's semi- major axis et^r, 
eccentricity ejv, and inclination i /y , with any speci- 
fied functiona l form . Following iMalhotral (|1993f ) and 
ILevison et alj (|2008l ) . we use the functional forms: 

e N = e exp(-f/T e ) 
i N = i exp(-t/Ti) 
aN = a f + (a - a/)exp(-i/r ) (1) 

where ao is the initial semi- major axis of Neptune, af — 
30 AU is the final semi-major axis, and T e , t;, and r a are 
the eccentricity damping timescale, inclination damping 
timescale, and migration timescale respectively. Our re- 
sults do not depend on the specific form of Eqn ([I]) . As 
we will demonstrate in Section [4j sometimes the instan- 
taneous rate of change of the variables (-, - , | ) is most 
relevant, while in other cases the total evolution matters 
most. We have verified these statements with integra- 
tions (not shown) using an alternative migration form 

- oc - oc I = constant. 

a e i 

3. PLANETESIMALS: OBSERVATIONAL CONSTRAINTS 
AND MODELED EVOLUTION 
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We model an initially unexcited disk of planetesimals 
that becomes today's cold classical population. In Sec- 
tion [3T] we present the observational constraints on the 
excitation of this population. In Section I3.2| we present 
an analytical model for the evolution of this planetesi- 
mal disk under the influence of Neptune, which we use 
to predict and interpret the results of numerical simula- 
tions. In Section 13.31 we justify directly modeling only 
Neptune instead of all four giant planets. 

3.1. Constraints from the observations of cold classical 

objects 

The cold classicals are a class of dynamically "cold" 
objects on low-eccentricity, low-inclination orbits, with 
positions starting at 42.5 AU, the region interior to 
which is unstable, and falling off quickly beyond 45 AU 
(|Kavelaars et alJ l2009h . We assume that today's cold 
classical KBOs are remnant planetesimals that formed 
in situ, and we use these terms interchangeably. Strong 
constraints can be placed on the dynamical history of the 
solar system by requiring that Neptune not disrupt these 
objects as it migrates outward on an inclined and/or ec- 
centric orbit. In this section, we discuss our specific cri- 
teria for "preserving" the cold classicals, which we will 
use to place constraints on the set of parameters (Section 
12.11) defining Neptune's orbital history. 

There is evidence that the classical KBOs have a bi- 
modal inclination distributio n ((Brown 200ll: lGulbis et all 
l2010HVolk fc Malhotrall2011l ). The cold classicals are de- 
fined as the class of objects with a distribution of inclina- 
tions i centered on a low inclination with a small width 
in inclination. The functional form of this distribution 
is ty pically modeled as a Gaussian multiplied by sini. 
One (jGulbis et al.l l2010) of the three proposed models for 
the de-biased inclin ation distribution differs substantially 
from the other two (|Brownll200ll : iVolk fc Maihotrall2011l ) 
with respect to the relative populations of the cold and 
hot classicals and the width of the hot (high i) compo- 
nent (Fig. [T] top panel). However, all three distributions 
are similar for the cold classicals (Fig. Q] bottom panel). 
The number of cold classicals per inclinati on bin falls off 
alm ost entirely by i = 6° for the models of lBrownl ()2001[ ) 
an d iGulbis et all (120101) a nd by % = 4° for the model 
of IVolk fc Malhotral (|201lD . Therefore we require that 
Neptune's dynamical history should not excite the cold 
classicals above an inclination of 6 degrees. 

The disruption criterion for the cold classical KBO ec- 
centricities is more subtle. The cold population is defined 
by its inclination distribution, not its eccentricity distri- 
bution. We could imagine a cold population of objects 
which have inclinations below six degrees but a uniform 
distribution of eccentricities. If this were true (and we 
will demonstrate that it is not), the initial planetesimal 
disk could be excited to arbitrarily large eccentricities. 
Moreover, the eccentricity distribution could be shaped 
entirely by the long-term stability of the KBO orbits. 
In this case, objects excited in eccentricity during Nep- 
tune's high-eccentricity period would be ejected from the 
system over billions of years. 

To test whether the eccentricities of the cold classi- 
cals are sculpted solely by stability, we compared ob- 
se rved cold classical object s to a sta bility map created 
bvlLvkawka fc Mukail (|2005D (Fig. [2]). iLvkawka fc Mukail 
(2005!) do not use proper elements, so we use the instan- 
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Fig. 1. — Thre e model s of th e un-biased in clination d istribution 
of the classicals: Brown (2001), dashed line; Gulbis et al.1 Q201Cf ). 
dotted line; and IVolk & Malhotra| II201H ), d ot-dashed line The 
biggest difference among them is that in the [Gulbis ct al] 112010 ) 
model, the hot (high inclination) component is a much less substan- 
tial portion of the total classical population and has a substantially 
smaller inclination width than in the other two models. The dis- 
tributions are very similar for the inclinations of the cold classicals 

feneous P orb 1 ital elements of the observed objects. We 
have confirmed that the features of the distributions we 
identify b elow are qualitatively the sam e using proper 
elements (|Dawson fc Murrav-Clavll2012l Appendix A.l). 
Because the hot and cold population overlap (Fig. [T]), we 
cannot definitively determine to which distribution any 
particular object belongs. However, for all three model 
distributions, less than 10% of objects with inclinations 
i < 2° are hot. We find that the eccentricities of i < 2° 
objects in the region from 42.5-45 AU are confined well 
below the survival limit. Consequently, we can conserva- 
tively constrain the dynamical history of Neptune: Nep- 
tune cannot excite the cold classical objects in this region 
above e = 0.1. 

Thus we can impose two conservative criteria for pre- 
serving the cold classicals: 

1. In the region from 42.5 to 47.5 AU, the inclinations 
of the cold classicals must not be excited above 

i < 6° 

2. In the region from 42.5 to 45 AU, the eccentricities 
must not be excited above e < 0.1. 

3.2. Secular evolution model for the planetesimals 

Excitation of the in situ planetesimal population by a 
perturbing planet on an inclined and/or eccentric orbit 
occurs through secular evolution. Here we present simple 
analytical expressions that we will use to predict and 
interpret the results of our integrations in Section 01 

3.2.1. Dynamics of secular evolution 

Due to forcing from NeptuneQ, both the eccentricity 
and inclination of a planetesimal undergo secular os- 
cillations on timescales of order a million years. The 

1 Other planets besides Neptune contribute to the secular forc- 
ing, but as a sim plification we consider only the effects of Neptune. 
See Section 13.31 for a detailed justification. 
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Fig. 2. — Plotted over the survival maps of Lvkawka & Mukai (2005) are the eccentricity (left) and inclination (right) distributions of the 
observed classical KBOs. The red squares are objects with i < 2° and are thus very likely cold classicals. The blue triangles have i > 6° 
degrees and are thus very likely hot classicals. The membership of any given purple circle (2° < i < 6°) is ambiguous. In eccentricity, the 
cold classicals (red squares) between 42.5-44 AU are confined to e < 0.05, well below the survival limit, while cold classicals between 44-45 
AU are confined to e < 0.1, also below the s urvival limit in thi s regio n. Classical objects are taken from the Minor Planet Center Database 
and classified by Gladman ct al. (2008) and Volk & Malhotra (2011). The yellow lines indicate the conservative criteria for preserving the 
cold classicals. 



planetesimal's total eccentricity is the vector sum of its 
forced eccentricity ef or cod, imparted by Neptune, and 
its free eccentricity ef rco , set by initial conditions (Fig. 
[3]). The efree vector precesses about the eforccd vector 
at the angular frequency <?kbo- The secular evolution 
of the planetesimal's inclination is analogous to - yet, 
to lowest order, separable from - the eccentricity evo- 
lution. The planetesimal's total inclination is the sum 
of its forced and free inclination, and the free inclina- 
tion precesses about the forced inclination. In the case 
of inclination, we can think of the particle's orbit being 
inclined by if ree with respect to the "forced" plane and 
processing about th e force d plane at the rate <7kbo- See 
iMurrav fc Dermottl (|2000f ) . Chapter 7, for a pedagogical 
presentation of secular evolution. 

The vector components of the planetesimal's eccentric- 
ity are h — e sin w and k — ecosro, where w is the 
planetesimal's longitude of periapse. Secular forcing by 
Neptune causes h and k to evolve as (to first order in e 
and cat): 



-0.2 
-0.2 




0.2 0.4 0.6 

k = e cosgj 



Fig. 3. — At a given time t, the total eccentricity (purple, solid) 
is the vector sum of ef orcc( j (blue, dotted) and ef ree (red, dashed). 
The vector ef rec precesses about ef orcet j at the rate </kbo> so a * 
time t, the ef ree vector has rotated by an angle (?kbo^ from its 
initial orientation with respect to ef orccd . 



h = e froc am(g KB ot + 0) + e fo rcod sin(njAr) 

k = e frcc COs(g K BC>£ + eforccd COs(n7Ar) 



(2) 
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where 



&a/a(«) 
a 



ejv, 



on 



hi 1 ) 

5KBO = 0:03/2 



(a) 



a 

m sun 4 



(3) 



The constants ef rco and ft are determined from the 
initial conditions. Here, vjn is the longitude of peri- 
apse of Neptune, ejv is the eccentricity of Neptune, and 
a is the ratio of Neptune's semi-major axis to that of 
the planetesimal, a, all of which are assumed to be con- 
stant. The functions b are standard Laplace coefficients. 
The secular frequency of the KBO is <?kbo, m-N is the 
mass of Neptune, m sun is the mass of the Sun, and 
n = (Grrisun/a 3 ) 1 / 2 is the planetesimal's mean motion. 

Similarly, the vector components of the planetesimal's 
inclination are p = isinfl and q = zcosfi, where is 
the planetesimal's longitude of ascending node. Secular 
forcing by Neptune causes p and q to evolve as (to first 
order in i and in): 



q = i free sin(-g K Boi 

p = i free COs(~5kBO* " 



-7) 
■7)- 



- ^forced Sin(fijv) 
forced COs(Q N ) 



where 



i N 



(4) 



(5) 



The constants Zf roc and 7 are determined from the initial 
conditions. Here, fijv is the longitude of ascending node 
of Neptune and ijv is the inclination of Neptune. While 
eforcod depends on a (Eqn. [3]), to first order if or cod is 
equal to ijv, regardless of the particle's semi-major axis. 

Three quantities are plotted in Fig. 0] as a function of 
a (for two different a N ): 1) ef orco d/eAr, 2) k OIce d/iN, and 



the secular period 



These values describe the am- 



plitude and timescale of the secular excitation of the in 
situ planetesimal population. For example, when Nep- 
tune is at 20 AU, a planetesimal at 45 AU has a forced 
eccentricity of 0.55gat and a forced inclination of iN- If 
the planetesimal begins with e = i = 0, it will reach its 
maximum eccentricity e = 2 x 0.55eAr = Lle^r and max- 
imum inclination i = 2i^ on a timescale of i „ 2 !7 =13 

Myr (i.e. half a secular oscillation period). 
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3.2.2. Secular excitation of the planetesimal disk 

The secular excitation of a planetesimal disk can be 
modeled using the expressions above. Here we consider 
that Neptune is temporarily on an inclined and/or ec- 
centric orbit after undergoing scattering on an effectively 
instantaneous timescale (i.e. much less than the secular 
timescale). Then Neptune imparts a forced eccentricity 
eforced and inclination iforced on the initially cold plan- 
etesimal disk. Thus each planetesimal has ef ree ~ ef orcec i- 
and if rcc ~ iforced- The planetesimal's total eccentricity 
and inclination, each a vector sum of the free and forced, 
now oscillates from the initial values of e(0) ~ and 
i(0) ~ 0, reaching maximum values of e = ef ree +ef orce d ~ 
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Fig. 4. — Parameters, to first order, governing secular evolution - 
using two different locations of Neptune, 20 AU (dashed line) and 
30 AU (solid line) - as a function of the planetesimal's semi-major 
axis. Top: Ratio of the planetesimal's forced eccentricity (dashed, 
solid) and forced inclination (thick, same for ajv = 20 AU, ajy = 30 
AU) to that of Neptune. To first order, the forced inclination does 
not depend on the position of the planetesimal relative to Neptune, 
while the forced eccentricity decreases with the planetesimal's semi- 
major axis. Bottom: Timescale of the secular evolution. 



2e forcc d and i = i flcc + i fo rcod ~ 2i force d, on a timescale 
set by the secular evolution rate <7kbo • Thus an initially 
cold planetesimal disk will become excited. 

3.3. The case for modeling only Neptune 's effects on the 
planetesimals 

We limit our parameter study to include only the 
planet Neptune. We have several reasons for choosing 
this approach. First, unlike previous approaches, we are 
not attempting to create a single model that reproduces 
the entire Kuiper belt in detail, but to constrain which 
histories of Neptune are consistent with a major qual- 
itative feature: the unexcited orbits of the cold classi- 
cals. Our constraints will feed into more detailed models. 
Second, restricting the parameter study to just Neptune 
drastically reduces the number of parameters, allowing 
us to thoroughl y explore the remaining parame ter space. 
Third, we find (|Dawson fc Murrav-Clavll2012[ ) that the 
primary effects of the other planets on the Kuiper belt 
are indirect: they alter the orbit of Neptune, which in 
turn affects the Kuiper belt. Our modifications to the 
Mercury 6.2 integrator (Appendix A) allow us to model 
any orbital evolution of Neptune without needing to in- 
clude the other planets. Fig. [5] provides one example 
pair of integrations, demonstrating that cold planetesi- 
mals evolve similarly in the presence of only Neptune and 
in the presence of all four giant planets. 

One main effect of other planets is to cause Neptune's 
longitude of periapse to precess. However, the disruption 
of the cold classicals is not significantly affected if Nep- 
tune's precession period is comparable to or longer than 
the secular excitation time (Fig. [5]). If Neptune pro- 
cesses very quickly, the cold classicals could be preserved 
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Fig. 5. — Snapshot at 0.5 Myr of test particles under the influence 
of Neptune alone (left) and all four giant planets (right). In both 
cases, the planetesimals begin at e = and Neptune has ajy = 
30, ejv = 0.2. Saturn, Jupiter, and Uranus are included in the 
integration on their current orbits (right) as an illustrative case. 
The black line is the same in both panels: the first-order predicted 
secular excitation of the planetesimals under the influence of just 
Neptune (Eqn. [2}, neglecting the effects of resonances. Outside 
of the resonances, the results are qualitatively similar. In this 
case (right), the forced precession period of Neptune due to the 
other giant planets (approximately 4.5 Myr) slow or comparable 
to the timescale of planetesimal secular evolution. Thus Neptune's 
precession has a negligible effect on the secular evolution of the 
planetesimals. 

outsid e our constraints, as proposed by iBatvgin et al.l 
(| 2 1 If) . In addition, strong interactions between Neptune 
and Uranus can cause Neptune's semi-major axis to oscil- 
late, leading to orbital chaos in the classical region, which 
can cause additional e xcitation. We explore these add i- 
tional complications in Da wson fe Murrav-Clavl (|2012f) . 

Here we do not explore the g$ and v% secular reso- 
nances, which over time remove low eccentricity and low 
inclination objects respectively. The location of these 
secular resonances will impose additional constraints on 
preserving the cold population. However, since they re- 
move unexcited objects, leaving excited objects behind, 
these resonances will not allow for parameters of Neptune 
that we rule out. 

3.4. Considerations for a massive planetesimal disk 

Throughout the paper, we treat the planetesimals as 
massless test particles. We parametrically model the or- 
bital evolution of Neptune caused by the planetesimals - 
Neptune's migration and the damping of its eccentricity 
and inclination - but do not explicitly consider the effects 
of planetesimal self-gravity. In this section, we consider 
some of these effects and how they might impact our 
conclusions. 

One important consideration is whether the transfer 
of angular momentum between Neptune and the disk is 
sufficient to excite a massive disk to the extent we as- 
sume for massless test particles. Since the total angular 
momentum of the system is conserved, the angular mo- 
mentum deficit (AMD) due to the secular excitation of 
the cold classicals cannot exceed Neptu ne's initial angu- 
lar momentum deficit (See[Laskar 1997, for a description 
of the concept o f AMD). Using the expansion of the AMD 
bv lHahnl (|2003l Eqn. 26a and 26b), the ratio of the total 
AMD of the cold objects to the AMD of Neptune is: 

AMp coid = ^coid / acoid ^ eg old + ig old ^ ^ 
AMD Ncptuno m N V a N e 2 N + i 2 N 

Since the factor y^jj^ is roughly unity, if the mass the 

cold classicals exceeds Neptune's mass, the AMD of Nep- 
tune would be too small to excite the cold classicals to 
ejy and ijv- Thus if the mass in cold classicals were large 
enough, the objects could remain at low eccentricities 



and inclinations even if Neptune's eccentricity and incli- 
nation were large. 

The outer solar system likely contained several tens 
of Earth masses at early times, comparable to Nep- 
tune's mass of 17M0. Such a massive disk is re- 
quire d to form large KBOs such as Pluto by coa gula- 
tion (jStern & Colwelll H997I: IKenvon fc Luul fl998h and 
to drive substant i al mi gration ((Fernandez fe Id 119841 : 
lHahn fc Malhotral I1999D and/or eccentricity and in- 
clination damping of N eptune's orbit. In contrast, 
iFuentes fc Holmanl (|2008T ) combined the results of several 
surveys to estimate that the mass of the current classical 
Kuiper belt is on ly 0.008 ±0.001M m , 40 % of which is the 
cold component (|Kavelaars et al.ll2009f ). 

We assume throughout the rest of this paper that the 
total mass - and thus AMD - of the cold classicals is 
smaller than that of Neptune (i.e. closer to today's 
mass). This choice implicitly assumes that the primor- 
dial disk was partially truncated or that the cold classi- 
cals were depleted before the era we treat in this paper. If 
the cold classical belt began with a low-mass - compared 
to the inner disk, responsible for Neptune's migration 
and for dynamical friction - truncation of the planetes- 
imal disk may explain why Neptune's migration halted 
at 30 AU ([Levison fc Morbidellil [2003D . Dynamical de- 
pletion would tend to excite the cold classical population 
to higher eccentricities than observed, but could be con- 
sistent if all excited objects were subsequently scattered 
out of the region from 42.5-45 AU. Thus this dynami- 
cal depletion would need to occur before the delivery of 
the observed hot classical population, which has higher 
eccentricities. During the era we consider in this paper, 
Neptune must deliver the hot objects without disrupting 
the cold ones. 

The constraints presented here will require revision if 
Neptune's primary sculpting of the Kuiper belt occurred 
while the disk was massive. For example, this could 
be true if the cold classical belt was depleted through 
collisions following Neptune's orbital evolution. Colli- 
sional depletion occurs through collisional grinding, fol- 
lowed by ejection of s mall grains by radiation forces, and 
IKenvon fc Luul (|1998f ) and lKenvon fc Bromlevl ([20011 ) ar- 
gue for this scenario in the context of coagulation models 
for the growth at KBOs. A break in the size distribution 
of the cold classicals at R = 25 — 50 km has b een at- 
tribu ted to collisiona l grind ing iPan fc Saril ([20051 ) . How- 
ever, [Nesvorn^eFiil] ([20111 ) argue that, given this inter- 
pretation, a large fraction of binaries with small compo- 
nents (R < 50 km) would have been disrupted. Such 
binaries are observed, implying that collisional grinding 
did not generate the break. Hence collisional grinding 
could only have substantially depleted the cold Kuiper 
belt if the majority of the belt's mass was initially se- 
questered in much smaller objects. 

Another potentially important effect is the propagation 
of sp iral density waves and torsion w aves in a massive 
disk (I Ward fc HahnlH998l : IBahnl [20031) . excited at secu- 
lar resonances and at the edge of the disk. In particular, 
spiral density waves tend to "smear out" the excitation 
of the eccentricities and inclinations of the planetesimals 
near secular resonances. Spiral density waves are poten- 
tially important in sculpting the dyna mical structure of 
the Kuiper belt - IWard fc Hahnl p98l) use them to place 
constraints on the mass of the Kuiper belt beyond 50 AU. 
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The exact behavior of the density waves depends on a va- 
riety of disk parameters - including the size distribution 
of planetesimals, the surface density profile in the disk, 
an d the t hickn ess of the disk - so we refer the reader 
to iHahnl (|2003[ ) for a detailed exploration. (For ease of 
comparison, note that assuming a surface density profile 
£ oc a^ 1 - 5 , the current mass of the cold Kuiper belt cor- 
responds to a £ roughly a factor of 3 smaller than the 
lowest surface density simulated by IHahnl (|2003l ).) With- 
out the influence of density waves, very large forced ec- 
centricities and inclinations (e — > oo, i — > oo) occur only 
close to secular resonances. The density waves spread 
the elevated eccentricity and inclination, confined to the 
secular resonance in the massless disk case, across the 
Kuiper Belt from 40-50 AU. However, our constraints do 
not include the extra excitation caused by secular reso- 
nances, which excite the cold classicals even more than 
we consider. Therefore, the constraints we will place on 
ruling out parameters of Neptune that cause excessive 
excitation will still hold and, as future work, the addi- 
tional excitation caused by secular resonances plus den- 
sity waves may place additional constraints. 

In the absence of secular resonances, waves launched 
at the disk edge could similarly smear eccentricities and 
inclinations. We assume that the eccentricities and in- 
clinations of planetesimals are not self-damped by spiral 
density waves launched at the disk edge. The eccen- 
tricities and inclinations of planetesimals could also be 
damped by dynamical friction caused by smaller, unob- 
served bo dies whose random vel ocities were damped by 
collisions ( Goldrei ch et al.l [2004) or by resonant relax- 
ation (|Tremainelll998f) . We assume that neither of these 
effects contribute significantly during the era of interest. 

Finally, the eccentricities and inclinations of the plan- 
etesimals may not be excited in the first place if the 
planctcsimal disk is massive enough to cause Neptune's 
eccentricity to quickly precess, lowering the forced ec- 
centricity of the cold classicals. In order for the disk to 
significantly affect the precession rate of the cold classi- 
cals, its mass must be at least that of Neptune. A sce- 
n ario involving a quick ly precessing Neptune is explored 
in lBatygin et al.1 (1201 lh and discussed in our companion 
paper. iDawson fc Murray- Clavl (|2012l ). 

4. RESULTS 

We present here the results of a set of integrations 
containing Neptune and a collection of test particles de- 
signed to represent the in situ planetesimals that become 
today's cold classical KBOs. The orbital evolution of 
Neptune is modeled as de scribed in Section 12.21 and the 
Appendix. In Section [4.1l we describe how we computa- 
tionally model the planetesimals. In Section l4~2l we char- 
acterize the secular excitation of the planetesimals under 
the influence of an inclined and eccentric Neptune and 
place limits on the eccentricity and inclination of Nep- 
tune. In Section 14.31 we determine the effect of damping 
of Neptune's eccentricity and inclination on the excita- 
tion of the planetesimals. We demonstrate that there are 
two regimes for damping Neptune's orbit - in the slow- 
damping regime, the planetesimals evolve to their initial 
free eccentricity and inclination, set by Neptune's ini- 
tial eccentricity and inclination, and in the fast-damping 
regime, their eccentricity and inclination are frozen at the 
value reached at the damping time. We also show that 



the evolution of the inclination and eccentricity can be 
treated separately. In Section l4~4l we consider the effects 
of Neptune's migration. Finally, we present two exam- 
ple scenarios in Section 14.51 demonstrating the principles 
we have derived: a pathological scenario in which the 
evolution of the Neptune's semi-major axis, eccentricity, 
and inclination occur each on a different timescale, and 
a realistic scenario consistent with preserving the cold 
population. 

4.1. Computational model of planetesimals 

We performed N-body integrations of an initially cold 
planetesimal disk under the influence of a dynamically- 
evolving Nept une using the Me rcury 6. 2 hybrid symplec- 
tic integrator (|Chamberslli999l ) with an accuracy param- 
eter of 10 -12 , a step size of 200 days, and user-defined 
forces and velocities imposing the migration and damp- 
ing of Neptune, a complete description of which is pro- 
vided in Section l2~2l and the Appendix. We modeled the 
initial population of planetesimals, which become today's 
cold classical KBOs, as 600 massless test particles with 
initial a evenly spaced between 40 to 60 AU and initial 
e = i = 0. We will use the terms "planetesimals" and 
"test particles" interchangeably from here forward. Al- 
though simulations that deliver KBOs into the classical 
region via scattering require ~ 10 4 planetesimals in each 
run, our immediate goal of understanding the conditions 
under which the cold classical population can be retained 
can be achieved with a smaller number, allowing us to 
explore a wider range of orbital conditions. 

Integrations probing solar system histories are typi- 
cally run for 4 Gyr in order to test the stability of the 
system under current solar system conditions. However, 
the dynamical histories of Neptune we test are not depen- 
dent on the long-term survival of the planetesimals, be- 
cause we have chosen observational criteria that are inde- 
pendent of the long-term evolution of the belt under the 
current configuration of the solar system (Secton lOT) . To 
save computation time, our integrations probe only the 
period in which Neptune undergoes planetesimal-driven 
migration and/or eccentricity damping and/or inclina- 
tion damping. 

4.2. Secular excitation of planetesimals under the 
influence of an inclined and eccentric Neptune 

In scenarios of the early solar system in which Nep- 
tune is scattered onto an eccentric and/or inclined orbit, 
Neptune can potentially disrupt an in situ population 
of planetesimals above the confined eccentricities and in- 
clinations we observe (Section 13. If) in the cold classical 
Kuiper belt today. Neptune forces the planetesimals to 
undergo secular evolution (Section l3.2|) . in which they os- 
cillate through high eccentricities and inclinations. The 
rate of secular evolution depends on the location of the 
planetesimal relative to Neptune (Fig. H]). Fig [5] shows 
snapshots of the secular evolution of test particles under 
the influence of Neptune when the planet is stationary at 
20 AU (top), migrates from 20 to 30 AU (middle), and 
is stationary at 30 AU (bottom). In their secular evo- 
lution, the planetesimals reach maxima e = 2ef or ced and 
i — 2if orcc d. In the middle panel, the evolution transi- 
tions from matching the top panel (early snapshots, left) 
to matching the bottom panel (late snapshots, right). 
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Thus, in the case of migration, we can estimate the sec- 
ular evolution using a constant on- When or in has 
not yet damped, we can estimate the secular evolution as 
if ajv were constant at the location where Neptune spent 
most time. Because we have chosen a migration law for 
which 2i£ decreases with time, planetesimals behave as 

if a at has always had the value ajy(t). We will discuss 
migration in detail in Section [4.41 

If Neptune's eccentricity and inclination are small (Fig. 
[7]) , the test particles remain below the observational limit 
of e < 0.1, i < 6° indefinitely. Since, to first order, 
iforccd is independent of a, the planetesimal's position 
relative to Neptune, the planetesimals will remain below 
i < 6° if i n < 3° . The forced eccentricity does depend on 
a (Fig. @J, decreasing with the planetesimal's distance 
from Neptune. Thus to satisfy our conservative criteria 
established in Section 13.11 which requires planetesimals 
with 42.5 < a < 45 AU to remain at e < 0.1, Neptune's 
eccentricity must stay below 0.09 at 20 AU or 0.06 at 
30 AU. However, if ejv and i^ damp due to dynamical 
friction with the planetesimals, the initial values of Nep- 
tune's eccentricity and inclination can be higher. In the 
next subsection, we discuss the effects of damping. 

4.3. Damping of Neptune's eccentricity and inclination 

The constraints we placed previously were on how high 
Neptune's eccentricity and inclination can remain indefi- 
nitely. During secular evolution, planetesimals reach ec- 
centricities up to 2ef orce( j and inclinations up to 2ef oxce d- 
However the planetesimals evolve to final values less than 
these maximum values, when Neptune's eccentricity and 
inclination damp. Here we refine the limits on Neptune's 
eccentricity and inclination in light of damping, consid- 
ering a range of timescales for the dynamical-friction- 
driven damping rates of ejv and ijy. The actual damping 
timescales depend on the surface density and size distri- 
bution of the planetesimals. 

4.3.1. Eccentricity and inclination can be treated separately 

To first order, the effects on the planetesimals of Nep- 
tune's inclination and eccentricity, including damping, 
can be treated independently. In first-order secular the- 
ory ( Section 13.21) . the evolution of a planetesimal's e is 
independent of Neptune's inclination i^ and of a plan- 
etesimal's i is independent of Neptune's eccentricity ejv- 
Thus we can place constraints on the evolution of Nep- 
tune's eccentricity without taking into account its incli- 
nation or vice versa. Fig. [8] illustrates the independence 
of the eccentricity and inclination parameters. 

4.3.2. Two regimes for damping: slow and fast 

The damping of Neptune's eccentricity ejv or inclina- 
tion in affects the secular excitation of the planetesimals 
in two regimes: slow and fast. We summarize constraints 
on Neptune's eccentricity and inclination in Table [TJ 

In the slow regime, e^v or i^ damps on a timescale T e 
or Tj, respectively, longer than the time (|27r/gKBo) f° r 
e or i to reach its maximum value. Because each plan- 
etesimal has an initial i(t = 0) = e(t = 0) = 0, its free 
inclination if ree and eccentricity ef rec are set by Neptune's 
initial e^(t — 0) and i^(t = 0). In this slow regime, the 
planetesimal's if rcc and ef ree are conserved. As e^r and ijq 
damp, the planetesimal's forced eccentricity ef orce d and 



TABLE 1 

Constraints on Neptune's eccentricity and inclination 



in 



no damping 



slow damping 



(a) 



1W 



fast damping 



-ejv < 0.1 



-ejv < 0.1 



,L ejvsin(g K BOTe) < 0.1 



2i N < 6° ijv < 6° ijv sin(g K BOT"i) < 6° 

inclination iforccd decrease, and the total eccentricity e 
and inclination i of the planetesimal approach ef ree and 

^free • 

Thus, in the slow regime, the planetesimals will evolve 
to a value below i < 6° if in < 6°. To satisfy our 
conservative criteria established in Section [3. 11 which re- 
quires planetesimals with 42.5 < a < 45 AU to remain 
at e < 0.1, Neptune's eccentricity must stay below 0.18 
at 20 AU and 0.12 at 30 AU. Thus in the slow damping 
case, there is a strong constraint (Table [TJ on the maxi- 
mum eccentricity and inclination at which Neptune can 
remain over long timescales. Note that these values are 
twice the values given in Section 1431 This is because in 
the slow damping regime, the planetesimals will evolve 
to ef rco = eforced (t = 0) and if rcc = iforccd(£ = 0), whereas 
if cn and in remain high indefinitely, the planetesimals' 
e and i continue to oscillate, reaching a maximum of 
2ef orcc d and 2i 

forced ■ 

If ejv and ijy damp in the fast regime, on a timescale 
shorter than the secular excitation time, if ree and ef ree are 
not conserved. The planetesimal's total e and i are frozen 
at the values they reach after approximately one damp- 
ing time. Fig. |H] illustrates the behavior of the planetes- 
imals in this regime in integrations in which Neptune's 
eccentricity or inclination damps. 

Neptune's orbit could have been even more eccentric 
and inclined than in the slow regime if eN and in damped 
quickly (Table [TJ right column) . If the planetesimal has 
not yet reached the maximum of its secular cycle after 
one damping time, instead of converging to the ef ree and 
ifr CC set by initial conditions, the planetesimal evolves 
to e fi na i = ef rflfi sin(t?KRor P .) and fe na j = if ree sm(#KBO n) 
(see lDawson fc Murray-Clavll2012l for a detailed expla- 
nation and justification). Consider again a planetesi- 
mal at 42.5 AU. If Neptune's eccentricity and inclination 
damp on a timescale of r = 0.32 Myr, the final eccen- 
tricity and inclination of the planetesimal are reduced 
by a factor of sin(gKBO T ) = 0.5 compared to the slow 
damping case. Thus the planetesimals will evolve to a 
value below i < 6° if i N < 6°/0.5 = 12°. To satisfy our 
conservative criteria established in Section [3. 11 which re- 
quires planetesimals with 42.5 < a < 45 AU to remain 
at e < 0.1, ejv must stay below 0.18/0.5 = 0.36 at 20 AU 
and 0.12/0.5 = 0.24 at 30 AU. 

We note that because "slow" and "fast" damping are 
defined relative to the secular evolution time, and be- 
cause the secular evolution time increases with the plan- 
etesimal's semi-major axis, it may be that the damping 
is "fast" for particles in the outer disk yet "slow" for 
particles in the inner disk. 

4.4. Migration of Neptune 

In the process of secular excitation of the cold classi- 
cals, migration alters a planetesimal's secular evolution 
timescale and forced eccentricity, which depend on a, the 



10 



Wolff et al. 



0.54 Myr 



2.73 Myr 



8.21 Myr 



27.3 Myr 



62.2 Myr 



186. Myr 




0.5 
0.4 

0.3 

0.2 

0.1 
0.0 

25 

20 

15 

10 

5 




0.5 
0.4 

0.3 

0.2 

0.1 
0.0 

25 

. 20 

1 

15 

' 10 

5 




45 50 55 
a (All) 



0.54 Myr 



45 50 55 
a (AU) 



2.73 Myr 



45 50 55 
a(AU) 



8.21 Myr 



45 50 55 
a (AU) 



27.3 Myr 



45 50 55 
a(AU) 



62.2 Myr 



45 50 55 
a(AU) 



186. Myr 



20 - 30 AU 
10 7 

0.2,i N =10 




0.54 Myr 

i N = 30 AU 
3n = 0.2J n =10 



2.73 Myr 



1 



8.21 Myr 



45 50 55 
a (AU) 

27.3 Myr 



li.'.'i 




45 50 55 
a(AU) 

62.2 Myr 



45 50 55 
a(AU) 

186. Myr 




:■_ : 



111 I I 1 \ 



ImIi,; 1 , i 



45 50 55 
a(AU) 



45 50 55 
a (AU) 



45 50 55 
a(AU) 



45 50 55 
a (AU) 



45 50 55 
a(AU) 



45 50 55 
a(AU) 



Fig. 6. — Snapshots of the orbital evolution of test particles beginning with e = i = from three different integrations (top, middle, 
bottom), each with ejv = 0.2 and ijv = 10° remaining constant. The top row of snapshots is from an integration in which Neptune's 
semi-major axis remains fixed at 20 AU; in the middle row, Neptune migrates from 20 to 30 AU on the timescale of r a = 10 Myr; and 
in the final row Neptune's semi-major axis remains fixed at 30 AU. The dashed line is the predicted first-order secular evolution of the 
planetesimals (Eqn. [2]and[4]|. The color of the planetesimals corresponds to the semi-major axis of Neptune, ranging from dark blue (20 
AU) to light blue (30 AU). The planetesimals in the integration that includes migration of Neptune (middle panel) matches the top panel 
(Neptune at 20 AU) at the beginning (left snapshot) and the bottom panel (Neptune at 30 AU) at the end (right snapshot). 



ratio of the planetesimal's semi-major axis to Neptune's 
(the forced inclination is independent a). For a plan- 
etesimal at 42.5 AU, the forced eccentricity is 40% larger 
when Neptune is at 30 AU vs. 20 AU and the secular evo- 
lution period 5 times shorter. Fig. [T0l shows the effect of 
migration on the secular evolution of the planetesimals. 
The excitation at a given time is well-modeled by the 
secular theory (Eqn. [5] and 2]) using Neptune at its snap- 
shot location. Recall that we are using a migration law 
that slows with time, so that Neptune spends increasing 
amounts of time at each, subsequent location. 

For retaining the cold classicals, what matters is the ra- 
tio of the migration timescale to the damping timescale. 
When this ratio is large ("slow" migration), Neptune ef- 



fectively damps at its initial position and we can model 
the secular evolution at this location (Fig. [TTJ top two 
rows). When this ratio is small ("fast" migration), Nep- 
tune effectively damps at its final position and we can 
model the secular evolution there (Fig. [TTJ bottom two 
rows). When the ratio is order unity, modeling the secu- 
lar evolution at the location Neptune reaches after half a 
damping time is a decent approximation (Fig. II H mid- 
dle row). However, we have not explored the regime in 
which r a r~j r e or r a ~ Tj in detail. Typically, we expect 
the damping and migration to occur in either the fast or 
slow regime, for reasons we will now state. The models 
in which Neptune is scattered from its loc ation of forma- 
tion to close to its current location (e.g. iLevison et al.l 
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Fig. 7. — A single snapshot, at 124 Myr, of test particles be- 
ginning with e = i = from two integrations (left and right) with 
the Neptune's ajv, ejv, ijv remaining constant. The integrations are 
identical except Neptune has ajv = 20, ejv = 0.09 in the left panel 
and ajv = 30, ejv = 0.06 in the right panel. Neptune has ijv = 3° 
in both panels. In both cases, the planetesimals are confine d to 
the low eccentricities and inclinations established in Section 13.11 
The dashed line is the predicted first-order secular evolution of the 
planetesimals (Eqn. [2] and Eqn. [1). 
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Fig. 8. — The eccentricity and inclination damping of Nep- 
tune have a separable effect on the planetesimals. Snapshots at 
22 Myr from three different integrations (top, middle, bottom) 
of test particles beginning with e = i = with Neptune at 
a N = 30, ejy = 0.2, iff = 10. In the top panel (blue), Neptune's 
eccentricity damps on a timescale of 0.3 Myr and its inclination 
remains constant. In the bottom panel (purple), Neptune's incli- 
nation damps on a timescale of 0.3 Myr and its eccentricity remains 
constant. In the middle panel (red), both the eccentricity and in- 
clination damp on a timescale of 0.3 Myr. In the middle planel, 
eccentricities of the planetesimals match the case in which just the 
eccentricity of Neptune damps (top) and the inclinations match 
the case in which just the inclination of Neptune damps (bottom). 



2008) can be considered fast migration. For extensive 
migration in a planetesimal disk, we expect the migra- 
tion timescale to be significantly longer than the damping 
timescale because the "random momentum" is a fraction 
of order of e or i the Keplerian momentum. (Though we 
note that Neptune only has to migrate some fraction of 
its semi-major axis.) 



4.5. Two example scenarios 

Consider two example dynamical histories for Neptune. 
In each, the planet starts at an initial location, eccentric- 
ity, and inclination and undergoes migration, eccentricity 
damping, and inclination damping to evolve to its cur- 
rent orbit. First we consider a pathological scenario in 
which the evolution of the Neptune's semi-major axis, 
eccentricity, and inclination occur each on a different 
timescale (Fig. [T2|) . Neptune begins at on = 20 AU, 
cn = 0.2, and in = 10°. On a timescale of r a = 3 
Myr, it migrates to its current location at 30 AU. The 
eccentricity damps on a longer timescale of r e =30 Myr 
and the inclination on a shorter timescale of Tj = 0.3 
Myr. For the inclination damping, the migration oc- 
curs in the slow regime (r a /rj > 1), and we model the 
damping as taking place at 20 AU. With Neptune at 20 
AU, the secular timescale of a planetesimal at 45 AU is 
27T = 24 Myr and the initial forced inclination of the 

9KBO J 

planetesimal is i forced = *jv = 10°. The value of incli- 
nation to which the planetesimal evolves is reduced by a 
factor of sm(gKBOTi) — 0.07. Thus the planetesimal at 
45 AU evolves to a final value of 0.7°. For the eccentric- 
ity, the migration occurs in the fast regime. Because the 
eccentricity damping timescale is much longer than the 
migration timescale (fast migration), we model the eccen- 
tricity damping as taking place at 30 AU. With Neptune 
at 30 AU, the secular timescale of a planetesimal at 45 
AU is 27r = 6 Myr and its initial forced eccentricity 

9KBO J J 

is eforccd = 0.78eAr = 0.16. The eccentricity damping 
is in the slow-damping regime, so the planetesimal's fi- 
nal eccentricity evolves to its initial forced eccentricity, 
e = 0.16. This predicted value and those for planetes- 
imals at other semi-major axes (orange curve, Fig. I12[) 
match the integrations well, validating our method. 

Finally, we consider a more realistic scenario consis- 
tent with preserving the cold population (Fig. HI?]) . Nep- 
tune begins at 26 AU, e» = 0.35, and %n — 14°. On 
a timescale of r a — 10 Myr, it migrates to its current 
location at 30 AU. Its eccentricity and inclination damp 
on the timescale of r e = t s ; = 0.3 Myr. In this slow 
migration regime, we model the damping taking place 
with Neptune at 26 AU. Thus the secular timescale of 
a planetesimal at 45 AU is 27r ~11 Myr, its initial 

L 9KBO J 

forced eccentricity is 0.7eAr = 0.24 and its initial forced 
inclination is ijv = 14°. The values of eccentricity and 
inclination to which the planetesimal evolves are reduced 
by a factor of sin^KBCe) = sin(<?KBOTi) = 0.16. Thus 
the planetesimal at 45 AU evolves to a final eccentric- 
ity of 0.04 and inclination of 2.4°, values satisfying the 
observational constraints established in Section 13.11 

5. CONCLUSIONS 

As a first step in a comprehensive study of the impact 
on the Kuiper belt of a wide range of possible dynamical 
histories of the outer solar system, we have performed 
a suite of numerical integrations probing the impact of 
the orbital evolution of a single a planet on a disk of 
planetesimals. We have presented the observational ev- 
idence for a population of dynamically cold objects in 
the Kuiper belt in the region from 42.5 to 45 AU that 
are confined to e < 0.1 and i < 6°. We argued that 
recent models of Kuiper belt sculpting - which explain 
many of the observed dynamical and physical proper- 
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Fig. 9. — Left: Snapshot at 44 Myr of the eccentricities and inclinations of test particles that began with e = i = 0. Neptune begins at 
a N = 30,ejv = 0, iff = 10. Then, in each of five different integrations, ipf damps on a different timescale: 10 Myr (purple, top), 3 Myr 
(blue, second from top), 1 Myr (green, second from bottom), and 0.3 Myr (red, bottom). When the damping timescale is longer than the 
secular evolution timescale (purple; blue particles interior to 47 AU), the planetesimals damp to their initial free inclination, which is set 
by Neptune's initial inclination. When the damping timescale is shorter (red; green; blue particles beyond 47 AU), the planetesimals are 
frozen at the inclinations they reached at approximately the inclination damping timescale. Right: Same for eccentricity damping, with 
ajv = 30, epf = 0.2, ijy = 0. The final eccentricity depends on the particle's location even in the slow damping regime (purple) because the 
initial forced eccentricity, the value to which the planetesimal's eccentricity converges, is a function of a, the planetesimal's semi-major 
axis relative to Neptune (Eqn. [3). In contrast, the forced inclination is independent of the particle's semi-major axis (Fig. [4{. Analytical 
curves are overplotted as orange dashed lines. 
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Fig. 10. — Snapshots of the orbital evolution of test particles 
(blue) beginning with e = i = from two different integrations 
(top, bottom), each with ejy = 0.2, ijv = 1-77° remaining constant 
and Neptune migrating from 20 AU to 30 AU. In the top model, 
the migration timescale is 10 Myr and in the bottom model the 
migration timescale is 1 Myr. The secular evolution model (dashed 
line) is computed at Neptune's location at the time of the snapshot, 
as if Neptune had been there during the entire secular evolution. 

ties of the Kuiper Belt - do not generate or preserve 
sufficiently low eccentricities for the cold classicals (e.g . 
iGome s 2003; Levison et aT1l2008HMorbidelli et al.ll2008l ). 



Our results have revealed several principles key to con- 
straining which orbital histories of Neptune are consis- 
tent with preserving the in situ planetesimal population 
at the low eccentricities and inclinations required by the 
observations: 

• If Neptune is scattered onto an eccentric and/or 
inclined orbit, it will secularly excite the eccen- 
tricities and inclinations of an in situ planetesimal 
population. 

• Planetesimals starting with e = i = reach eccen- 
tricities and inclinations up to twice their forced 
eccentricity and inclination, on timescales given by 
Eqn. ([21) and (g]). 

• As Neptune's eccentricity and inclination damp, 
the planetesimals evolve to their final eccentricities 
and inclinations. If the damping timescales r e and 
Ti of Neptune's orbit are slow compared to the secu- 
lar evolution time, a planetesimal evolves to its ini- 
tial free eccentricity and inclination, which are set 
by Neptune's initial eccentricity and inclination. If 
the damping is fast compared to the secular evo- 
lution time, the planetesimal's e and i effectively 
freeze at the values they reach after one damp- 
ing time, reaching final values of ef ree sin(<7KBO' r e) 
and eforcod sin(gKBO T 'i) respectively. See Table[T]for 
constraints on Neptune's eccentricity and inclina- 
tion in the two damping regimes. 
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Fig. 11. — Snapshots at 54.7 Myr from 10 different integrations. Left: Effect of the migration timescale when Neptune's inclination 
damps. In each case, Neptune has ejv = 0, an initial iff = 10, and an inclination damping timescale of r, = 1 Myr. In this fast-damping 
regime, Neptune's inclination damping timescale is shorter than the planetesimals' secular evolution timescale, so the inclinations of the 
planetcsimals are frozen at approximately the values they reach at Tj = 1 Myr. In the top row, Neptune's semi-major axis stays fixed at 20 
AU (infinitely slow migration). In rows 2-4, Neptune's migration timescale is 10 Myr, 3 Myr, and 1 Myr respectively. In row 5, Neptune's 
semi-major axis stays fixed at 30 AU (infinitely fast migration). The color of the planetesimals corresponds to Neptune's semi- major axis 
at the time of rj/2 = 0.5 Myr, half the damping timescale, ranging from dark blue (20 AU) to light blue (30 AU). Because Neptune's 
semi-major axis sets the secular evolution time, the inclination that a given planetesimal reaches before it is frozen depends on the migration 
timescale relative to the inclination damping timescale. In rows 1-2, the model (dashed line) is computed with apf = 20 AU. In the middle 
row, it is computed with ajv = 24 AU. In the rows 4-5, it is computed with ajv = 30 AU. Right: Same for eccentricity damping: in each 
case Neptune has ejv = 0.2, an initial ijv = 0, and an eccentricity damping timescale of r e = 1 Myr. 



• The effects of Neptune's: 1) eccentricity evolution, 
and 2) inclination evolution, on the planetesimals 
can be treated separately to first order. 

• At a given location in the planetesimal disk, the 
secular excitation timescales and forced eccentric- 
ity (but not the forced inclination) depend on Nep- 
tune's location, which is altered by Neptune's mi- 
gration. When Neptune's migration is slow rel- 
ative to the damping time, the secular evolution 



effectively takes place at Neptune's initial loca- 
tion. When Neptune's migration is fast relative to 
the damping time, the secular evolution effectively 
takes place at Neptune's final location. 

From these pr inciples, it is evident that the three mod- 
els described in iLevison et al.l ((2008) w °uld not be able 
to retain a cold classical population. Neptune begins 
with an eccentricity of 0.3 and a semi-major axis of 27.5 
AU (run A and C) or 28.9 AU (run B) and damps on 
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Fig. 12. — In this example integration, the semi-major axis, inclination, and eccentricity of Neptune evolve on three different timescales: 
T a = 3 Myr, Tj = 0.3 Myr, r e = 30 Myr. The initial conditions of Neptune are ajv = 20 AU, ejv = 0.2, ijy = 10°. In each snapshot, the 
color of the planetesimals corresponds to Neptune's semi-major axis, from 20 AU (dark blue) to 30 AU (light blue). For the eccentricity, 
the model (dashed line) is calculated at the location of Neptune in the particular snapshot. For the inclination, because the inclination 
damping is fast compared to the migration, the model (dashed line) is calculated at Neptune's initial location of 20 AU. 



a timescale of 1 Myr (run A and B) or 3 Myr (run C), 
in the slow migration regime. When Neptune is at 27.5 
AU, a planetesimal at 42.5 AU has a forced eccentricity of 
e = 0.76ejv = 0.23 and a secular evolution timescale of 6 
Myr. The final eccentricity of the planetesimal evolves to 
is reduced by a factor of sin(gKBO r e) = 0.85 for a damp- 
ing timescale of 1 Myr (fast damping) and is not reduced 
for a damping timescale of 3 Myr (slow damping). Thus 
the final eccentricity of the planetesimal is 0.19 (run A) 
or 0.23 (run C), well above the observational limit. When 
Neptune is at 28.9 AU (run B), a planetesimal at 42.5 
AU has a forced eccentricity of e = 0.79eAr = 0.24 and 
a secular evolution timescale of 5 Myr. The final eccen- 
tricity the planetesimal evolves to is reduced by a factor 
of sin(gKBO T e) = 0.97 for a damping timescale of 1 Myr. 
Thus the final eccentricity of the planetesimal is 0.23, 
well above the observational limit. According to the con- 
straints established in our paper, any of these three initial 
conditions for Neptune could retain the cold classicals if 
Neptune's eccentricity were to damp more quickly. 

Based on these principles, we can place robust con- 
straints on Neptune's dynamical history. For example, 
in the regime of slow damping and slow migration, if 
Neptune's initial semi-major axis is 20 AU, then its ec- 
centricity must stay below 0.18 and inclination below 6 
degrees. If Neptune's initial semi-major axis is 30 AU - 
or if it migrates quickly, relative to the damping time, 
from its initial location to 30 AU - then its eccentricity 
must stay below 0.12. In the case of fast damping - on a 
timescale shorter than a planetesimal's secular excitation 



time - the initial eccentricity and inclination of Neptune 
can be even higher. 

Having established these principles, we complete 
our parameter study i n tw o companion papers 
iDawson fc Murray-Clayl (j2012) and Dawson and 
Murray-Clay (2012b), in prep. In the first companion 
paper, we consider more generally the constraints on 
Neptune's dynamical history from the eccentricity dis- 
tribution of the classical KBOs and incorporate several 
other important effects in the constraints from the cold 
classicals: 

• A more accurate model for secular excitation that 
includes higher-order terms. 

• The greatly increased secular frequency near 
Neptune's mean-motion resonances, which places 
strong constraints on Neptune's semi-major axis 
when its eccentricity is high. 

• The effects of the other giant planets, including pre- 
cession of Neptune, whic h can reduce a plan etesi- 
mal's forced eccentricity (|Batygin et al.ll20TT| ) , and 
oscillations in Neptune's semi-major axis, which 
can create a chaotic sea in the Kuiper belt region. 

Then we combine these constraints for retaining the cold 
classicals with constraints for creating the hot classical 
population and identify which regions of parameter space 
of Neptune's dynamical history can produce the hot clas- 
sical population without disrupting the cold population. 
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Fig. 13. — Example integration for a region of the parameter space of Neptune's orbital evolution that is consistent with producing 
the observed cold classical population. Neptune migrates from 26 to 30 AU on a timescale of r a = 10 Myr. Its initial eccentricity and 
inclination are 0.35 and 14° respectively and both damp on a timescale of r e = t; = 0.3 Myr. In each snapshot, the color of the planetesimals 
corresponds to Neptune's semi-major axis, from 26 AU (medium blue) to 30 AU (light blue). The model (dashed line) is calculated at the 
location of Neptune in the particular snapshot. 



In the second companion paper, we place constraints on 
Neptune's inclination and inclination damping time. 

Our parameterization of Neptune's orbital history al- 
lows us to constrain which orbital histories are consis- 
tent with maintaining the cold classical population. By 
combining the observational constraints (Sec. I3.1[) with 
the principles derived in this paper, we can immediately 
check whether a particular set of initial conditions for 
Neptune's semi-major axis, eccentricity, and inclination, 
and the rates of its migration and eccentricity and in- 
clination damping, are consistent with maintaining the 
observed dynamically cold population, without perform- 
ing computationally expensive integrations. A picture, 
both qualitative and quantitative, is emerging of what 
dynamical histories of Neptune allow a promising general 
scenario - Neptune's delivery of the hot classicals from 
the inner disk to classical region, where the cold popula- 
tion has formed in situ - to be consistent with observed 
low eccentricities and inclinations of the cold classical 
population. Barring fast precession of Neptune's orbit, 



which we do not consider in this work, the existence of 
the cold classical population implies that Neptune could 
have spent only a limited time at high eccentricity and/or 
inclination during its dynamical history. 
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APPENDIX 

MODIFICATIONS TO THE EQUATIONS OF MOTION FOR THE EVOLUTION OF NEPTUNE'S ORBITAL ELEMENTS 

In this appendix, we follow ILee &: Peald ((2002) to derive modifications to the equations of motion of an orbiting 
body to allow any form of evolution of the body's orbital elements. We use these modified equations to model the 
early solar system evolution of Neptune's orbit caused by interactions with the other giant planets and with the 
planetesimal disk. Directly modeling Neptune, the other giant planets, and tens of thousands of massive planetesimals 
would be computationally prohibitive. Instead, we model Neptune alone, with its orbit evolving as it would under the 
influence of the other planets and the planetesimal disk, including migration, damping of its eccentricity, and damping 
of its inclination. However, the equations we derive can be used to implement any type of orbital evolution, not just 
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migration and damping. For example, in a companion paper (|Dawson fc Murrav-Clavl[20Tl) . we use the modifications 
to cause Neptune's semi-major a xis to oscillate, a s it would due to resonant interactions with Uranus. 

We modify the Mercury 6.2 (|Chambersl [1999) N-body integration code to allow for arbitrary orbital evolution 
by adding extra terms to the equations of motion at each timestep. We add these terms via a user- defined force, 
alread y supported by Mercury 6.2, and an additional velocity modification described below. We note that lLee fc Pealel 
(2002)^ use their derived modifications for a non-symplectic code and use a symplectic implementation when employing 
a symplectic algorithm. However, we do not add the extra terms in a symplectic manner, even for a symplectic 
integration algorithm; thus these modifications must represent small perturbations on the planet. We determined that 
a timestep of 200 days was sufficiently small for our models, yielding the same results as smaller step sizes. 

We make these modifications as follows: 

1. At each timestep, Mercury 6.2 calls the user-defined force routine. We have modified this routine to return not 
only a user-defined acceleration but a user-defined velocity. We also modify the routine to only apply these 
corrections for Neptune, not for the planetesimals. 

2. The user-defined force routine calculates the user-defined veloc ity addition al terms from Eqn. (|A9|) - (|A1 1|) below 
and the user-defined acceleration additional terms from Eqn. (fA12|l - ((AT4]) below. 

3. The acceleration additional terms are used to advance the velocity of Neptune and the velocity additional terms 
to advance the position of Neptune, in addition to the usual gravitational forces. 

To derive these additional terms, we follow iLee fc P eale (2002) but instead of holding the orbital inclination i 
constant, we allo w it to vary with time 

Equation (|A1|) gives the position vector components (x, y, and z ) as a function of the orbital elements (compare 
to Eqn. A5 of ILee fe Peale). The variable r is the distance of the planet from central body, fl is the longitude of the 
ascending node in the xy plane, oj is the argument of periapse, / is the true anomaly, and i is the inclination. 



x = r cos il cos (oj + /) — r cos i sin £1 sin (oj + f) 
y = r sin £1 cos (oj + f) + r cos i cos f2 sin (oj + f) 
z = r sin i sin (oj + /) 



(Al) 



The velocity vector components are given in equation (|A2[) (compare to Eqn. A6 of ILee fc Peale! ) . A dot over a variable 
indicates its derivative with respect to time. 



x = cos £1 [f cos (oj + f) — rf sin (oj + /)] — sin £1 [r cos i sin (oj + f) + rf cos i cos (oj + f) — zi] 
y = sin £1 [r cos (oj + /) — rf sin (oj + /)] + cos £1 [f cos i sin (oj + /) + rf cos i cos (oj + /) — zi] 
z = r sin i sin (oj + /) + rf sin i cos (oj + /)] + ri cos i sin (oj + f) 



(A2) 



Below are the x-components of the "velocity" and "acceleration" additiona l terms used t o update the body's position 
and velocity, respectively, at each timestep (compare to Eqn. A3 and A4 of ILee k. Pealef ). The variable a is the semi- 
major axis and e is the orbital eccentricity. 
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Similar expressions can be derived for the other coordinates. It should be noted that these are simply the additional 
terms to the equations of motion resulting from the orbital evolution - d, e, and/or i - and do not describe the 
overall motion of the system. All coordinates, velocities, and accelerations must be calculated with respect to the 
central body. In order to compute these partial derivatives, we must first define the variables r, r and rf following 
IMurrav fc Dermottl (|2000f) : 



q(l-e 2 ) 
1 + e cos / 



VT^e 



-.e sin f 



rf = 



VT 



= (1 + ecos/) 



(A5) 



Ne xt we calculat e the partial derives of r, f, and rf with respect to the orbital elements a, e, and i (compare to Eqn. 
A7 of ILee fc Pealel ): 
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(A6) 
(A7) 

(A8) 



Combini ng equations (IA6I) through (| A8|) with equation (|A3|) yields an expression for change in position (compare to 
Eqn. A8 of lLee fc PealeT) : 
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We can now use these additional terms to update the position of a body undergoing any arbitrary orbital evolution in 
a, e, and/or i, including migration, eccentricity damping, and inclination damping. Below we combine (|A6I) through 
(IA8I) with eq uation (|A4j) to obtain the acceleration terms used to update the body's velocity (compare to Eqn. A9 of 
ILee fc Pealef ): 
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— cos (u + f) — sin (uj + /) 



cos f2 
COS ft 



2o. 

dr d(r f ) • dr z 

— cos i sin (u> + f) H — cos i cos (lu + f) — i— — 

de de de r 



z H — z 
i 



— sin i sin (to + f) 
2a 

dr 

— sin i sin (d + f) - 
de 



— rf r • 

- — — sin i cos (w + /) H — i cos i sin (w + /) 
2ft a 

d(rf) dr ■ 

— - — sin i cos (w + /) + — i cos i sin (w + /) 
ae ae 



r cos i sin (u) + f) + r/ cos i cos (w + f) + r- cos i sin (u + f) — iz 

i 



(A12) 



(A13) 



(A14) 



Arbitrary functions for the evolution of th e or b ital e lements can be assigned to the a/ a, e/e and i/i terms. With 
the additional terms derived above (Eqn. IA9I - |A"14)) . the full motion of a body undergoing arbitrary evolution in 
semi-major axis, eccentricity, and/or inclination can be implemented. When the eccentricity e reaches a small value, 
machine round-off error can cause the eccentricity damping terms to damp the eccentricity below 0, leading to failure of 
the integrator. Therefore did not apply the eccentricity damping terms if e < 10 _4 .If the orbit of the body precesses, u> 
is no longer constant (apsidal precession) and/or f2 is no longer constant (nodal precession). Since in this paper we do 
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not consider precession (see lBatygin et al.ll201lt iDawson fc Murrav-Clavl F2012'). we leave the additional modifications 
to the equations of motion to account for precession for future work. 
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